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Different approaches to solve the spinor-spinor Bethe-Salpeter (BS) equation in Euclidean 
space are considered. It is argued that the complete set of Dirac matrices is the most appro- 
priate basis to define the partial amplitudes and to solve numerically the resulting system of 
equations with realistic interaction kernels. Other representations can be obtained by per- 
forming proper unitary transformations. A generalization of the iteration method for finding 



o 
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the energy spectrum of the BS equation is discussed and examples of concrete calculations 



are presented. Comparison of relativistic calculations with available experimental data and 
with corresponding non relativistic results together with an analysis of the role of Lorentz 



H , boost effects and relativistic corrections are presented. 

A novel method related to the use of hyperspherical harmonics is considered for a repre- 
sentation of the vertex functions suitable for numerical calculations. 



I. INTRODUCTION 

The interpretation of many modern experiments requires often a covariant description of 
two-body systems. This is due either to the high precision that calls for an inclusion of all 
possible corrections to a standard (possibly non relativistic) approach, or due to the high energies 
and momenta involved in the investigated processes. In the subatomic field the most obvious 
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examples are the properties and structure of the deuteron and, to some extent, the mesons, 
if the these are treated as quark antiquark systems. Also the charm degree of freedom is a 
recent topic of analysis in heavy ion experiments. The CBM experiment of the future FAIR 
project at GSI will investigate highly compressed dense matter in nuclear collisions with a beam 
energy range between 10 and 40 GeV/u. An important part of the hadron physics project is 
devoted to extend the SIS/GSI program for the search of in-medium modifications of hadrons 
to the heavy quark sector providing a first insight into charm-nucleus interaction. Thus, the 
possible modifications of the properties of open and hidden charm mesons in a hot and dense 
environments is matter of recent studies. 

Within a local quantum field theory the starting point of a relativistic covariant description 
of bound states of two particles is the Bethe-Salpeter (BS) equation. Due to known difficulties 
related to the analytical properties of the BS amplitude in Minkowski space, the procedure 
of solving the equation turns out to be a rather involved task. Up to now, the BS equation 
including realistic interaction kernels has been solved either in Euclidean space within the lad- 
der approximation (see e.g. Refs. hjl- 5fl and references quoted therein) or adopting additional 
approbations of the BS Ration itself ^.The procedure of solving „„ m e ri cah y the BS 
equation has been revisited in Ref. [131 ] . and a reduction of the BS equation to an equation of 
the Light Front form has been proposed in Ref. [9(. Moreover, a detailed investigation of the 
solution and the properties of two spinor particle in a bound state within Light Front Dynamics 
have been reported in Refs. [a, Il4[ ]- 

The present paper represents a brief review of mathematical and numerical procedures of 
solving the BS equation in Euclidian space. Advantages and disadvantages of different methods 
depending of the feature of the attacked problem are discussed. Section ILTl surveys basic features 
in solving the BS equation. 

In Subsection III Al it is demonstrated that the Dirac representation of the BS amplitude 
with subsequent use of the iteration method for solving the system of integral equations for the 
partial amplitudes, is the simplest and most appropriate way to obtain a numerical solution. 
However, the physical interpretation of results of calculations within such a representation is 
not so transparent. The spin angular harmonics basis, as discussed in Subsection III Bl provides 
a more transparent interpretation of the partial amplitudes in terms of familiar notations used 
in non relativistic approaches. Obviously, in practice, one can combine the two representations; 
the numerical solution may be obtained within the Dirac representation and then, by a unitary 
transformation, transferred to the spin angular basis. A covariant form, i.e. a frame independent 



representation of the BS equation is presented in Subsection III CI The possibility to obtain 
approximate analytical expressions of matrix elements within the BS formalism is discussed in 



Subsection III Dl in the context of the one-iteration approximation procedure. A generalization 
of the iteration procedure to find the energy spectrum of the BS equation is presented in Section 
IIIII The use of the hyper spherical harmonics basis to interpret the BS amplitude in form of 
one-dimensional numerical arrays, suitable for solving the BS and Dyson-Schwinger equations 
for quark-antiquark mesons, is discussed in Subsection IIVI 

II. SOLVING THE BS EQUATION 

The BS amplitude for two spinor particles A and B is defined as [15[ 

*(xi, x 2 )o/3 = (0\T^(x 1 )^(x 2 )\AB) = e- iPX ^(x) aP , (1) 

where \AB) denotes the bound or scattering state of the AB system, tp (a?i) and ip B (x2) are 
the fermion field operators in the Heisenberg representation with spinor indices a, /3, and x = 
x\ — X2 is the relative coordinate. The amplitude ([T]) obeys the BS equation which in the ladder 
approximations reads 

^ v- 2 f H S(pi) 76 (l)\W)7b(2)5(P2) , 9 s 

^(P) = iy.9b / TTTXi 7 75 T-- ' ( 2 ) 

k J ( 27 (p-py-nl + ie 

where b denotes the type of exchanged mesons (e.g. n, p, u>, a, 5 and r\ for the deuteron problem); 
76 are the corresponding nucleon meson vertices; the nucleon momenta are px,2 = P/% ± P, P, 
and S(p) is the free nucleon propagator. The amplitude (pQ) is redefined as ^ = — $ Uc, where 
U~c is a charge conjugation-like matrix (Uc = ^7270) or Uc = 7173) S(q) = Uc S T (q) Uc- 

From Eqs. ((TJ) and (|2j) it can be seen that in Minkowski space the BS amplitude is a rather 



complicate object, being a 4 x 4 matrix in the spin space and having cuts [la] and poles along 
the real axis of the relative energy in configuration space. Usually, instead of solving directly the 
matrix equation (|2|), one decomposes the BS amplitude \I/ over a complete set of basis matrices 
and solves the equation for the coefficient functions of such a decomposition. Generally, these 
coefficient functions form a system of four-dimensional integral equations. In order to reduce 
their dimension the coefficient functions are also decomposed over a complete set of orthogonal 
functions in configuration space. Such a procedure is mathematically correct if the original 
amplitude is an analytical function of its arguments, which is not the case in Minkowski space. 
To obtain a completely regular equation one performs the Wick rotation [16] to Euclidian space 
where the unknown quantities (i.e. the partial BS amplitudes) are defined along the imaginary 



axis of the relative energy. The Mandelstam technique [17( relates the matrix elements of 
observables in physical processes with the BS amplitudes in Euclidian space. However, often, 
the BS amplitudes are needed directly in Minkowski space, so that an algorithm for analytical 
continuation of the Euclidian solution in the whole plane of the relative energy po or a numerical 
procedure of inverse Wick rotation back to Minkowski space are demanded. These procedures 



are not yet well established and are still under consideration by many theoretical groups [35(. In 
the present paper, we focus our attention on mathematical and numerical methods of solving the 
BS equation in Euclidian space with realistic interaction kernels in the ladder approximation. 
The problem of finding the solution in Minkowski space is not considered here. 

A. Dirac representation 

The most appropriate choice of the basis for the representation of the amplitude \l/(p) in 
spinor space is the complete set of Dirac matrices 7, (i = 1 . . . 16), which contains a minimal 
number of 4x4 matrices and, consequently, significantly simplifies calculations of corresponding 
matrix elements and traces. In this representation the amplitude becomes 

tf (p) = lil> a (p) + tbVpCp) + 7a.C(p) + Tfe7/*C(P) + <vO), (3) 

where the subscripts of the coefficient functions ipi refer to the transformation properties, i.e., 
these coefficients can be scalar (s), pseudo-scalar (ps), vector (v), axial- vector (a) and tensor (t) 
functions, respectively. Further, in order to eliminate the angular dependence, the coefficients 
i\)i are decomposed correspondingly over the spherical Ylm and vector angular Y^ M harmonics. 
For instance, in the 3 Si — 3 D\ channel (that is for the deuteron case) the amplitude reads 

®d(p) = T7 J75^plYiAf + 757oV>alYlM - 1pvl(j, Y\ M ) - 1pa075h,YlM) 

IpI l 

-^275(7, Y? M ) - 2#?i7o(7, Y 1m) + #to(7, [7, Y° M ]) + # t2 ( 7 , [7, Y? M ])} • (4) 
By placing eq. (HJ) in to eq. ([2]), calculating the corresponding traces and carrying out angular 



integrations one obtains a system of eight two-dimensional integral equations of the form 
ippi(jPo,\p\) = Ki 
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^aiOo, |p|) = Ki { -2m N poij pl + 



ml } -p> + 2pl-- A M 2 d 



Ci- (6) 




- 2p |p| I ~7=>l>a0 ~ \l^a2 | } , (7) 

^i(Po,|p|) = Ki| }• (8) 

where Ki is an integral operator whose explicit expression depends on the type of the exchanged 
meson. For instance, for scalar mesons it reads 

Ki^(P0, IpI) = -i&Dfa, |p|) J dP '°^^ 2 Q L (y)^ L ( p ' 0l |p'|), (9) 

and 

* P) = (S^F« (10) 

is the scalar part of two propagators. In Eq. Q, Ql{v) are the Legendre functions with the 
argument y = — — ■ — U ^/ i - — — with \lb as the mass of the exchanged meson. Analogous 
expressions can be obtained for exchanges of pseudo-scalar and vector mesons. 

The resulting system of integral equations ([5])- (|5]) , after performing the Wick rotation, is 
solved by using the iteration procedure. Note that, since the eigenvalue M^ (i.e. the deuteron 
mass) enters nonlinear ly in Eqs. in ©-((HI) (cf. Eq. (110j) ). the iteration procedure is highly 
prevented. Usually, one fixes the binding energy at the experimentally known value and searches 
the solution relative to the coupling constants g which enter linearly in Eqs. ([5|)-([H|) and for which 
the iteration procedure can be employed. 

We solved by this method the BS equation for the deuteron with a realistic kernel within the 
one-boson exchange interaction with six exchange mesons |J|. The numerical solution consists of 
eight partial amplitudes in the form of two-dimensional arrays indexed by |p| and po, which can 
be further used to analyze the properties of the deuteron 1 191 and to compute matrix elements of 
different processes [19|]-[20(]. It should be noted that, besides advantages (simple Clifford algebra 
with Dirac matrices, fast convergence of the iteration procedure) within the Dirac representation, 
the physical interpretation of the partial amplitudes (J3J)- (SJ) and matrix elements of observables 
in terms of familiar nonrelativistic quantities, Lorentz boost effects and comparisons with other 
approaches (light front dynamics [a, |l2(, Gross equation [6|] etc.) are difficult. 



B. Spin angular basis 

The basis of spin-angular harmonics T^- (p) [l|, y, [2l| allows for a more transparent interpre- 
tation of the solution. It is constructed from the complete set of solution of Dirac equation for 
free nucleons 

ndPWc = (-) P1+P2 i L E (Lm L Sms\JM) ( l -p 1 l -p 2 \Sm s )Y LmL {p) Ufc{p) Uff{-p){ll) 

where p = p/\p\, and the spinors U^\{pj) are the solutions of the free Dirac equation for the 
particle "i" with spin projection p and positive or negative p-spin [211 ]; a denotes other quantum 
numbers of the system; a = {LSJp±p2}, where L is the relative orbital momentum, S the 
total spin and J the total momentum of the system. Often, for a the spectroscopic notation 
a = 2S + 1 L^ lp2 is employed [21] . The explicit expressions for the spin-angular harmonics can 
be found, e.g., in Ref. |l9j]. Evidently, the Dirac and the spin-angular harmonics basis are 
connected via a unitary transformation, presented explicitly in Ref. [19]. Now the strategy can 
be formulated as follows: the BS solution and matrix elements of observables are obtained within 
the Dirac representation, then by using the unitary transformation to the spin-angular basis, 
results are rewritten in a more familiar form in terms of spectroscopical partial amplitudes. This 
provides a more clear physical meaning of the obtained expressions and allows for a detailed 
comparison with other relativistic and non relativistic approaches. Moreover, with partial BS 
amplitude within the spectroscopical classification one can "a priori" estimate the contributions 
of different terms. So, those amplitudes which have a direct analogue in the non relativistic 
approach shall provide the main contribution at low energies. Other partial waves can contribute 
only at higher energies and higher momentum transfer. For instance, for the deuteron the eight 
components can be classified in three groups: 1) the main components with two positive values 
of the p spin, ^ 3>g ++, (j> 3D++y which are the relativistic generalization of then familiar S and D- 
waves of the deuteron 2) four "negative" P- waves, with one positive and another negative p-spin: 
4> i p+- iip-+, <^_a p+- , 4> 3 P -+ , and 3) two negligible small components with both p-spins negative 
4>3s — ' ( i >3 D — [19]- Usually, the later ones are disregarded in most calculations. As an example 
of the obtained solution we present in FigQ] the positive components S ++ and D ++ of the BS 
vertex (solid lines) and compare them with the non relativistic results obtained from Schrodinger 
equation with Bonn and Paris potentials (dashed and dotted lines, respectively) and with results 
provided by the relativistic three-dimensional Gross equation [6]. As expected, at low values 
of the relative momenta p all approaches provide basically the same results. Differences occur 
at large values of p, where the validity of non relativistic approaches becomes questionable. 
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FIG. 1: Comparison of the deuteron S (left panel) and D-waves (right panel) obtained within differ- 
ent approaches: Bethe Salpeter |4j (solid lines), Gross equation [6j (dot-dahes lines), non relativistic 
Schrodinger approach with Bonn (dashed lines) and Paris (dotted lines) potentials. 

With this solution we have analyzed different reactions with the deuteron in electromagnetic 

B 



20, 



22j, weak 



23] and hadronic processes 



24 



261 ] with the result that calculations within the 



BS formalism provide a much better agreement with experimental data in comparison with non 
relativistic approaches. 

C. Covariant representation 



It should be noted that the non relativistic Schrodinger approach belongs to the so-called 
"instant-form" of dynamics [271 ] . i.e. the solution of the Schrodinger equation is obtained in 
the center of mass system of the two-body system; being invariant only relative to the Lorentz 
transformation which leave the plane t = const unchanged, it can not be boosted to another 
system of reference. This means that in non relativistic calculations the Lorentz boost effects in 
the wave function of a moving particle are always ignored. Probably this can be justified only at 
low momentum transfer (cf. e.g. Fig. [1]). Contrarily, the BS amplitude, Eq. (JTJ), is a covariant 
object by definition. Consequently, by applying the Lorentz transformation to the amplitude 
found in the rest system of the two-body system, one can obtain the BS solution in any system 
of reference. However, mathematically such a procedure is quite involved (see, e.g. Ref. |19l]). 
and in real calculations it becomes rather difficult to separate and investigate pure boost and 
other relativistic effects. In order to avoid such difficulties it is tempting to find a covariant 



representation of the amplitude, i.e., a representation within which the partial amplitudes are 
frame independent. Such a representation can be accomplished if one observes that, for instance 
in the deuteron case, one can form a new covariant basis of 4 x 4 matrices from the available 



four independent matrix structures, £, (p£,)I, S(pi 



\-i 



P P 

+ p — m, and S{p2) _1 = — — p — m, 



2 ' yi ' 2 

where £ is the deuteron polarization four-vector. Then the BS amplitude can be presented in 



the form 



*d(P,p) 
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in 



h 7 £ + h t 



(pO 



in 



m 

P/2 — p + m 



in 



(12) 



where the hi are the sought covariant partial amplitudes. Obviously, these new amplitudes can 
be related to the non covariant ones by projecting (|12p onto the corresponding basis. Explicit 
relations of covariant amplitudes hi in terms of spin-angular components can be found in Ref. 
[19l | . Below we present an illustration of how the covariant amplitudes (fT2|) can be used to 
analyze the process of backward elastic pD scattering. This process has been investigated in the 
non relativistic Schrodinger approach within which the wave function in initial and final states 
has been taken for the deuteron at rest. The corresponding non relativistic helicity amplitudes 



and the cross section can be written as 24| 

2 
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&NR 



w{P lah ) (2V2u(P lab ) - w(P lab )) P h 



Cnr = (u(P lab ) + ^^) (<Plab) ~ V2w(P lab )) Pi 



labi i 
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MPi 
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P, 



labi 



r/r ' V/,! -Siu^P^ + w^P^fp, 
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(13) 

/„,„ (14) 



V2 v \/2 )~ la01 dn 

where Pi ab is the momentum of the nucleon in the laboratory system, u and w are the two 
components of the deuteron wave function and the partial helicity amplitudes are defined as 

T = A d M e MI ) + B (n£ M )(n(+ ,) +iC(a- [£ M x ^+ ,]) + W (an)(n ■ [£ M x ^+ ,]). (15) 

In Eq. (fT5|) £ M is the polarization 3-vector of the deuteron at rest, T is the total amplitude 
of the process. This amplitude can be found by calculating the corresponding one-nucleon- 
exchange Feynman diagram and using the covariant BS amplitudes (|12p which allow to write 
explicitly the dependence on the polarization vector £ M and, consequently, to determine the 
partial amplitudes A,B,C and T>. The result of a calculation, with preserving only the main 



" + +" components, is [24] 

dvo n firAfrt , , , T . 2/ri ^2 ^ A , P* , 29P* SU»; 



dn 3 mPlab) + *^)) 2 ^( 1 + ^ + ^ + ^ + ...), (16) 

A = i? a6 (MA*) - ^|^) 2 -C(P^), (17) 

Bo = PL \^D{PlaJb) (2^2^ s{Plab) + ^D(Plab)) C(P lab ), (18) 



Co = PL (ysiPlab)- * ^ ) (^ S (Plab) + V2^ D (P lab )) C(P lab ) 



(19) 



V = - PL ^fn{Piab) (*s(Plab) ~ * D ^ afe) ) iCCJ^), (20) 

where ^s^d are the main, " + +", components of the deuteron (see Fig. [[]), and C{Pi ab ) = 
I 1 + x 1 ^ + 5 lo i +•••)• It can be seen that the analytical expressions for relativistic and non 
relativistic helicity amplitudes are rather similar, except for the factor £(Pi ab ) which depends 
on the velocity of the deuteron and obviously describes the Lorentz boost effects. Other pure 
relativistic corrections appear in the cross section from the contribution of the negative P- waves 
[24| . In Fig. [21 we exhibit results of calculations of the differential cross section with emphasis 
on different contributions. The notation in Fig. [2] is quite self-explanatory. One can conclude 
that at large values of the initial energy the relativistic effects increase and become significant 
at Pab > 0.35 - 0.4 GeV/c. 

D. The one-iteration approximation 

As mentioned above, the basis of Dirac matrices turns out to be the most appropriate repre- 
sentation in solving numerically the BS equation. Transitions to other bases can be accomplished 
by performing the corresponding unitary transformations. In such a way one can obtain numer- 
ical solutions of the BS equation in representations, which are most relevant for the considered 
problem and one can investigate numerically different aspects of the process. However besides 
numerical investigations, it is useful to have a mathematical method for analytical estimates 
of the considered effects. Remind, that often in non relativistic approaches one takes into ac- 
count relativistic corrections by calculating additional diagrams with meson exchange currents, 
nucleon-antinucleon pair currents etc. It is intuitively clear that some of these corrections must 
be included already in the relativistic amplitude via negative P waves and Lorentz boost effects 
which are absent in the non relativistic calculations. An estimate of the role of P waves can be 
accomplished by solving the BS equation in the so-called one-iteration approximation (OIA). 
The spirit of the method follows from the observation that if, when solving numerically the BS 
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FIG. 2: The spin averaged cross section da/dCl for the backward elastic dp scattering in the center of 
mass as a function of the momentum of the detected proton in the laboratory system. Left panel: Dashed 
line - contribution of the positive BS waves; long-dashed line - Lorentz boost effects, solid line - full BS 
calculations; dotted line-results of the non relativistic approach with the Bonn potential wave function. 
Experimental data are from [281 ]. Right panel: dashed line - contribution of the one-nucleon exchange 



mechanism only; dotted line - contribution of the triangle diagram (A-isobars); solid lines - the full BS 
calculations with P waves (a) and without P waves (b). 

equation by iterations, the trial function is close to the exact solution, the latter one can be 
obtained by performing only one or two iterations. 

As illustrated in Fig. [H at low values of the relative momentum p the exact solution coincides, 
to a good accuracy, with the non relativistic solution of the Schroedinger equation. This is a hint 
that, using the non relativistic deuteron wave function as trial function, the number of iterations 
in solving the BS equation can be strongly reduced. Moreover, at low and intermediate values 
of p the first iteration with the non relativistic solution as trial function can serve as a good 
approximation of the exact solution. 

As an example, let us consider the BS equation for the deuteron partial vertices in the 
spin-angular harmonics basis (fTT|) : 



G a (k , |fc|) = ~ X> 2 / ^y ° W a p(k,p) ^(po, \p\), 



(21) 



Wafi(k, P ) = J2 



I in 



Qi(z) 

27r|p||fc| 



dn k d%Y lm {p) YCnik) Tr[r+(-fc)76r /3 (p) 7b ] 



The OIA method implies: i) in (|21|) all "negative" waves are put equal to zero, ii) the fco and po 
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dependencies in the partial kernels W a p(k,p) and in vertices G 3s ++ and G 3£) ++ are disregard, i.e. 
W a p(k,p) ~ W a/ g(|fc|, |p|), G ++ (po, |p|) ~ G ++ (0, |p|), iii) the trial function is taken as solution 
of the non relativistic Schroedinger equation for the deuteron, and iv) the negative (P) waves 
are obtained only by one iteration in Eq. (|2ip . For instance, for the pseudo-scalar isovector 
exchange we get 



W p -+_> s++ 



-W p+ - 



>s++ 



W r _ 



>D++ 



-w p +^ D++ 



V2W p+ ~ 
= y/2W p + 
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V2W P -+ 
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>D++ = \/2Wp i -+_ >D+H 



■NV2 
= N 



Qo(y)\k\-Qi(y)\p\\(?2) 
-Qi(y)\p\ + Q2(y)\k\]m) 



where J\f = — y/3/(2\p\\k\Ek), Ev. = Vk 2 + m 2 ~ m, m being the nucleon mass. In coordinate 
space one has 

p 2 dp fif*,^ t N l _ 1/n , {\ iTr n -j X f j ^fs(f) 
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2\p\\k\ 

p 2 dp 

2\pW\ 



\k\Q (z)-\p\Qi(z) *s(\p\) 



\0\Qi(z)-\k\Q 2 (z) * D (\p\) 
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-e-^il + nrMkr), 



(24) 
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^d{t) 



-/jjt 



(l+pr^kr), (25) 



where fi stands for the pion mass. Then the result for the negative P-waves for the pseudo-scalar 
one-boson exchange reads 



* P+ - 



(n + ~) 



(Plab) 



2 2^/3 



-fir 



dr (1 + fj,r) hirPlab) [N u u(r) + N w w(r)] 



(26) 



where u{r) and w(r) are the non-relativistic deuteron wave functions in the coordinate repre- 
sentation, and g 2 ~ 14.5 is the pion-nucleon coupling constant. The normalization factors are 
N u = V2 (1) and N w = -1 (v2) for P 3 {P± ) waves. With this negative-energy waves one 
may estimate the origin of the relativistic corrections computed within the non-relativistic limit 
as additional contribution to the impulse approximation diagrams, such as meson exchange 
currents and NN pair production currents. As an example we present here the relativistic 
corrections to the partial amplitude A (|15p 
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(27) 



which is similar to the expressions obtained in non-relativistic evaluations of the so-called "catas- 



trophic" and pair production diagrams in electro-disintegration processes of the deuteron 



29] 



and also similar to the resu 
the elastic pD processes 



24 



ts of a computation of the triangle diagrams usually considered in 



30, 



311 ] . Another example concerns the calculation of the effects 



of final state interaction in reactions of deuteron break-up with two nucleons in the final state 
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FIG. 3: The cross section for the reaction pD —>■ (pp)n. Left panel: dotted and dashed lines represent the 
calculations within the nonrelativistic and relativistic impulse approximation, respectively; the two solid 
lines are the results which include the final state interaction in the (pp) - pair in a : 5o state calculated 
within the one-iteration approximation. The lower and upper solid lines have been obtained with and 
without cut-off form factors in the BS equation. Right panel: the relative contribution of the deuteron 
relativistic S (dotted line) and D (dashed line) components to the total cross section. Experimental data 
are from 32 1 



with low relative energy [32]]. Exact calculations require also the solution of the BS equation 
for the scattering state in the continuum. At the considered energies one can avoid the involved 
procedure of solving the BS equation in the continuum by employing the OIA in the same 
manner as described above. In Fig. [3] we present results of calculations of the cross section for 
the reaction pD — > {pp)n within the OIA for the BS amplitudes and a comparison with non 
relativistic results and available experimental data. It can be seen that relativistic calculations 
provide a much better agreement with data. 



III. EXHAUSTING METHOD 



Notice that the above mentioned methods are higly suitable for solving the BS equation by 
iterations, which provide the solution of ground state of the system. This is quite acceptable 
in most calculations for the deuteron and light mesons. However, heavier mesons (e.g. J/\£, 
D mesons) possess an entire mass spectrums with various, experimentally known, transitions 
to the ground state. Consequently, one needs a generalization of the iteration procedure which 
would allow to find the energy spectrum of the excited states and the corresponding amplitudes. 



The "exhausting" (depletion) method 33] is the one within which the ground and excited states 
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of the system can be found by solving the BS equation by iteration. The main idea of the 
method can be illustrated in the case of the scalar BS equation with scalar interaction. The 
corresponding equation reads 

/d A p' 1 
7^4 7 A2 TT~ VJMip'), (28) 

(27TJ 4 [p-py - fi z + te 

where the scalar propagator is S(p) = l/(p 2 —m 2 + ie) and JM are the total angular momentum 
and its projection. It is more convenient the introduce the BS equation for the vertex function, 
¥Jm(p) = S{p\) Gjm(p) S{p2)- Then in Euclidian space one has 

Wp) = a JH piyj m + pl _ lMl? l plMl]{(f _„ f + ^ <W>, (29) 

V v ' 

K(p,p';M B ) 
where A = g 2 . Obviously, at fixed Mb Eq. (|29p is a Fredholm type eigenvalue equation relative 



to Xu = gjj with eigenfunctions G JM (p)- The set of eigenvalues Ajj = Ajj(M^) is also called 
the spectrum of (|29p relative to the coupling constant. It is straightforward to establish a corre- 
spondence between the mass spectrum of the equation and its spectrum on coupling constant; 
it is sufficient to find the spectrum Ajj(Mb) at all allowed values of Mb, < Mb < 2m and 
to inverse the problem, i.e. to fix A at its actual value and to find the corresponding value of 
Mb- In such a way, each Aj receives its particular value of Mb, which forms the mass spectrum 
of the equation. From the rigorous Hilbert-Schmidt theory of symmetric integral equations it is 
known that: i) the spectrum of Eq. ([29]) is discrete and real with a nondecreasing sequence of 
Kj |Ai| < IA2I < ••• < |A n | . . ., ii) the corresponding eigenfunctions Gjm(p) an d the correspond- 
ing amplitudes tp l j M {p) form a bi-orthogonal basis f d 4 p (p' L J M {p) Q 3 ji M i{p) = SijSjj'Smm' within 
which the kernel K(p,p'; Mb) can be represented as 

00 _. 

K(p,p>; M B ) = Y,Y,— g JM(p) ^M(P')- (30) 

i=l JM lJ 
For the scalar equation the partial amplitudes have a particularly simple form, Gjm(p) = 
t4tGj(p)Y j m{^p) , where for brevity the notation p = (po, \p\) has been introduced. Then 
the partial vertices G l j{p) read 

G'AP) = >»j [tip' k j (p;p';M b ) g;0), (31) 

kJ{p ^ Mb> = w? [(3 +ri -iU)'+ri*a QAv) - (32) 
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Notice that the partial vertices Gj{p) together with the partial amplitudes (p j also form a 
bi-orthogonal basis, hence the partial kernels K (ft;p'\ Mb) can also be written in the form 

oo 1 

K J (p;p>; M B ) = Y, T ~ &M ifW). (33) 

Now we are in the position to formulate the exhausting method for the Eq. (|29p . Consider the 
problem of finding the lowest eigenvalue \\j at fixed J and Mb- Schematically, Eq. (|3ip read 

-L &j = K J G 3 . (34) 

Let now choose a trial function x(p) and form a sequence of functions if %, (if ) 2 x, •••, (if ) n X- 

oo 

The trial function % can be decomposed on the complete set of functions Gj, x = Yl c * Gj, 
where q = J" dp <fj{p) x(P) an d from Eq. (|34"|) one gets 

oo 1 OO -j oo _. 

^ J X = E C ^^ G ^ (if J ) 2 X = £> ^ G^, • • • , (K J Yx = Y,Ci^G*j.{^) 
£i A * J i= i A u 7=t A ^ 

At each n-th iteration one considers the ratio {K J ) n x/{K J ) n ~ 1 x, which at large enough n 

provides the lowest eigenvalue of the spectrum: 

on 



(K J Yx 


Z^Ci jw- G l j 

i=l u 


( K 3)n-^x 


oo 

Yl C % yn-\ G j 






i=2 x 7 n— >-oo 



(36) 






i=l »J i=2 

Equations (|35p and (j36|) demonstrate how one finds the ground state solution, i.e. the 
eigenvalue Au and the eigenfunction G\j can be found. Now, to find the next excited eigenvalue 
A2J, one defines a new kernel if" 7 ' 1 for which the new BS equation has the same spectrum Xu 
and eigenfunctions Gij (i = 2...n) except the first state i = 1. For such a kernel, A2J and G23 
would serve as ground state solution. The new kernel which obeys these conditions is 

00 

K J ^(p,p'; M B ) = K J (p,p>; M B ) - — Gj(p) <fitf) = £ -j- G W) <fiW)- (37) 

Au ^ Xu 

It can be seen that all eigenfunctions Gj of the kernel K , except Gj, are the eigenfunctions of 
the new kernel K J)l with the same eigenvalue spectrum Ajj. Obviously, the lowest eigenvalue 
of K J ' 1 is A2J, which can be found by the same iteration procedure, i.e. via Eqs. ([35]) and ([36]) . 
Analogously one finds the next solutions i = 3, 4 . . . from the recurrent kernel 



1 °° 1 

K J > n (p,p>; M B ) = K^&p'; M B ) - — G n ; (p) tftf) = £ — GS(p) ¥>S(pO- (38) 

A " J i=^i AiJ 
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FIG. 4: The coupling constant g vs. the value of the mass Mb of the bound system for two values 
of the total angular momentum J = 0, 2. The solid curves correspond to the ground state, while the 
dashed lines are for the first excited state. Intersections of the horizontal line g — 10 with the respective 
curve determine the energy spectrum of the system: Mb — 1-75 GeV , Mb — 1.94 GeV for J = and 
M B ~ 1.97 GeV, M B ^ 1.98 GeV for J = 2. 

As an example of application of the method, results of calculations of the eigenvalues g = j 
as functions of the binding mass Mb are presented in Fig. U] for first two ground states of a 
system of two equal masses particles with total angular momentum J = and J = 2 (solid lines) 
and two excited states (dashed lines), respectively. It is clear that, if the coupling constant is 
known from independent experiments, the mass spectrum of the system can be easily recovered 
(see illustration in Fig. 2] for g = 10). 



A generalization of the method for the spinor-spinor BS equation can be found in 331 ] . 



IV. HYPERSPHERICAL HARMONICS 



The considered methods are convenient for solving numerically the BS equation by an iter- 
ation procedure. The obtained solution presents a set of two-dimensional arrays of the partial 
amplitudes ipi, cf. Eq. Q. In practice, such a representation of the solution can cause difficul- 
ties in specific calculations which require two-dimensional interpolations of such arrays, e.g. in 
attempts to determine analytical parametrizations of the numerical solution, in performing the 
inverse Wick rotation, etc. This is particularly awkward to use such solutions if one solves the BS 
equation for qq bound systems (mesons) in nonperturbative QCD. In the later case, the BS equa- 
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tion must be solved simultaneously with the Dyson-Schwinger equation to generate dynamically 
the masses of the constituent quarks. Numerical solutions in form of one-dimensional arrays 
would be much more appropriate for tthese problems. This can be achieved by decomposing 
the partial amplitudes over the three-dimensional hyper spherical harmonics basis in Euclidian 
space 



/92m ( n _i_ 1 \( n _ Al/12 
Wx,M)=*m(x)*WM); Xm(x) = d - l ( n ^ + 1) ; ^ l xC l n + _)(cosx), (39) 



IPI • A - . Ll 



where C^jJ are the usual Gegenbauer polynomials, cosx = ^~r', sinx = ^r) p = \l Pq + P 2 - 
The BS vertices (amplitudes) are decomposed over the complete set of spin-angular harmonics. 

G(po,p) = ^g a (po,\p\)T a ( P ) (40) 

a 

with the coefficients g a (po, \p\) as 

00 

9a(P0,P) = ^2 9 "(P} Z jlm{XpiQpAp)- (41) 

j,n,m 

Note that, in order to be able to carry out the v p , 9„, <p p - integrations analytically, the spin- 
orbital basis must be slightly redefined (see Ref. 34| for details). The scalar part of the 
interaction kernel is also represented in terms of hyperspherical harmonics 

1 



Op - k) 2 + n 2 

ram 



271 " 2 5Z ~^TJ V n(P> k)Znlrn{Xp, Op, 4>p) KlmiXk , Ok, (t>k) 



4 /A, -A- 



By inserting Eqs. ()40|) -f)42f) into the BS equation one obtains a system of one-dimensional 
integral equations for the partial amplitudes 



9\{po,\p\) = ^2gi(p)x 2j -2,o(xp); 92{po,\p\) = 5Z 5f 2(p)^2j-2,o(xp); (43) 

00 00 

9a(po,\p\) = Y^9i(p) x 2j-i,i(Xp); 94(pq,\p\) = ^gl(P) x 2j,i(x P ), (44) 

where 
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/dkk 6 
87r2(2j _ 1} V 2j -2(P, fc)EE A )i "(*) C(k), (45) 

r\ 71^—1. 7Yl^ _L 

/,f L3 4 oo 

|^ F 2i _!(p, k) E E 4J?(*) »»(*), ( 46 ) 

n ^ n=lm=l 

flj(p) = V&4 / 2f , ^(p, fc) E E tt C(*)-, (47) 

Q W "I" J n = lm=1 

where 6j and A!?^ 1 are known coefficients, whose explicit expressions can be found in Ref. [34]. 
It can be seen from Eqs. (|45p-(|47p that the partial amplitudes within the hyperspherical basis 
are functions of only one variable. Instead, the obtained system of integral equations is infinite. 
At first glance, this may cause unresolvable problems. However, in practice, when solving the 
system (|45|) -(|45 |) a good convergence can be achieved by considering only the few first terms 



in (|43p and (]4"4"|) , usually up to j = 3 — 4 [34| . In Fig. [5] we present an example of the use of 
hyperspherical basis in solving the BS equation for the 1 Sq channel. The obtained solution can 
be fitted by a simple formula 



9i(p) - 



P 2 



+ 6? 



V 



3 1 4 n « 2 '- 2 

E T^^T ( 4 8) 



which is extremely useful in further analyses of the solution in Minkowski space. The fitted 
parameters a,j and bj can be found in Ref. [34]]. 

V. SUMMARY 

In this mini-review we briefly considered different methods to solve the BS equation in Eu- 
clidian space. It is shown that the complete set of Dirac matrices is the most convenient way 
to decompose the BS amplitude and to find the partial amplitudes by iteration procedure. The 
transition to other representation, e.g. to the spin-angular harmonics basis or to the covariant 
representation can be accomplished by performing the corresponding unitary transformation. 
Examples of numerical solution for the deuteron and applications for the calculation of matrix 
elements of observables are presented, which illustrate the role of relativistic corrections and 
Lorentz boost effects. 

A generalization of the iteration method to find the energy spectrum of the BS equation is also 
discussed. It is shown that the Hilbert-Schmidt bi-orthogonal basis can define a mathematical 
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FIG. 5: Functions g{,j — 1, . . . , 4, eq. (|4"5|) - (l4"4"|) . Full squares correspond to g\\ full circles to g\ multiplied 
by 10; triangles to g\ multiplied by 100; open circles to g\ multiplied by 1000; the solid lines correspond 
to the fitted functions g{ by formula (|48j) . The overall normalization constant is arbitrary. 



method, known as the "exhausting method", to find numerically the energy spectrum and the 
corresponding partial amplitudes for the excited states of two-body systems. 

The hyperspherical harmonics basis is shown to be the most appropriate one in finding the BS 
solution in form of one-dimensional numerical arrays. This is extremely convenient in analytical 
continuation of the solution back to Minkowski space and in solving simultaneously the Dyson- 
Schwinger and BS equation for the qq bound states (mesons) in nonperturbative QCD. 
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